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Directional pulse propagation in beam, rod, pipe, and disk geometries 



(N 

o 

(N 
-i— > 
O 

O 

in 

(N 



C/3 . 

o- 

• <— I . 
^— > 

q. 

o 

Of 



> 

OV 

O' 
(N 



X 



Blackett Laboratory, Imperial College London. 



P. Kinslerf] 

Prince Consort Road, 

th 



London SW7 2AZ, United Kingdom. 



(Dated: Friday 26 th October, 2012) 

I derive directional wave equations useful for pulses propagating in beam, rod, pipe, and disk 
geometries by using a cylindrical coordinate system; the scheme works equally well for either long 
multi-cycle or single-cycle ultrashort pulses. This is achieved by means of a factorization procedure 
that conveniently generates exact bi-directional and first order wave equations after the selection 
of propagation direction - either axial, radial, or even angular. I then discuss how to reduce these 
to a uni-directional form, and discuss the necessary approximation, which is essentially a paraxial 
approximation as appropriately generalized to the specific geometry. 



I. INTRODUCTION 



Directional decompositions of wave equations |l|-|5j arc 
a powerful method of developing pulse propagation equa- 
tions valid down to the ultrashort and few cycle regime. 
Done correctly, these provide exact bi-directional forms, 
which can be systematically approximated in a "slow evo- 
lution" limit into a uni-direction form of great practical 
use. Usually, the emphasis is on the role of nonlinearity 
or dispersion, and on modelling the propagation of beams 
or pulses in a linear geometry. 

A feature of the decomposition is the choice of a refer- 
ence evolution, containing as much of the detail of the 
exact propagation as possible, with the rest left as a 
(hopefully) perturbative "residual" term that couples the 
forward and backward waves. Fortunately, even if the 
residual is not that weak, the backward contribution is 
exceedingly poorly phase matched, which greatly extends 
the validity of the approximation. Nonlinearity is typi- 
cally left as a residual term, as is diffraction. 

Here my intention is to take an alternative path, and 
provide the basics of bi- and uni-directional propagation 
models for non-Cartesian geometries, although at first I 
restrict myself to the cylindrical case relevant to beams, 
rods, and disks; others can be easily generated as re- 
quired. In particular, although the result of factorization 
for directed axial propagation will likely look familiar, the 
radial and angular cases are more interesting. In particu- 
lar, such propagation models have potential applications 
for disk, ring, and other whispering galley optical res- 
onators [Tc|] . My focus on geometry rather than disper- 
sion or nonlinearity, means that it is the diffraction terms 
that are of more interest than nonlinearity and disper- 
sion; for these the "slow evolution" criteria amounts to a 
generalised paraxial approximation. Although the overt 
focus is on optical pulse propagation, since the Helmholtz 
wave equation used as a starting point is useful in many 
fields (e.g. for acoustic pressure waves), the results here 
have potential for wider application. 



II. THEORY 

Most optical pulse problems consider a uniform and 
source free dielectric medium. In such cases a good start- 
ing point is the second order wave equation, which results 
from the substitution of the V xH Maxwell's equation 
into the Vx£ one in the source- free case (see e.g. p|). 
Further, assuming linearly polarized pulses, we can use a 
scalar form. Defining V 2 = d 2 + d 2 + d 2 and d a = d/da, 
we can write the wave equation as 



1 



— d t e * fi-k 



E(t) = Q. 



(1) 



Here I have suppressed the spatial coordinates for nota- 
tional simplicity; in fact we have E(t) = E(t, r) and the 
total polarization Q(E,t,r); also r = (x,y,z). Note that 
a full expansion of the various possible components of Q 
is given in [0 , but in summary it can contain nonlinear- 
ity, dispersion, and free current effects - and potentially 
even magnetic nonlinearity. It can even allow for non- 
Helmholtz behaviour, such as that present in some acous- 
tic wave models || ; and if adapted can generate tempo- 
rally propagated wave equations instead of the spatially 
propagated ones described here |g. 

Of the potential complications, we here include the 
isotropic linear material response terms in a reference 
wavevector fc 2 (w) = e(uj)fi(u!)u! 2 . Thus, in the frequency 
domain, we can write 



-k 2 (uj)] E{t) = -Q. 



(2) 



In most descriptions of pulse propagation we will want 
to chose a specific propagation direction and then denote 
the orthogonal components as transverse behaviour. Of- 
ten this process uses Cartesian x,y,z coordinates (see 
e.g. H, but here I show how directional techniques can 
be applied in alternative geometries. 

I now factorize the wave equation, a process which, 
while used in optics for some time Jfl2j has only recently 
been used to its full potential [|J ||, • Given a wave 
equation of the form 
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K 2 \ E = — Q, 



(3) 
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with Q now also including the non-<9 2 derivative terms, 
we can see that the LHS of eqn. (||) is a simple sum of 
squares which might be factorized, indeed this is what 
was done in a somewhat ad hoc fashion by Blow and 
Wood in 1989 ([3. Since the factors are just d z =F iK, 
we can see that each (by itself) would generate a for- 
ward directed wave equation, and the other a backward 
one. Leaving basic mathematical detail to the appendix, 
a rigorous factorization procedure JjJ [l3| allows us to de- 
fine a pair of counter-propagating Greens functions, and 
so divide the second order wave equation into a pair of 
coupled counter-propagating first order ones. 

Counter-propagating wave equations suggest counter 
propagating fields, so I split the electric field up accord- 
ingly into forward (E + ) and backward (E~) parts, with 
E = E + + E~ . The coupled first order wave equations 
are 



d z E ± = ±zKE ± ±^|. 



(4) 



The RHS now falls into two parts, which I term the un- 
derlying and residual parts |15[ . First, there is the iKE^ 
term that, by itself, will describe plane-wave like propa- 
gation in the simplest cases. Second, the remaining part 
oc Q which can be called "residual" terms. These residual 
contributions, here containing the transverse derivatives 
®x + ®y> account for the discrepancy between the true 
propagation and the underlying propagation. Although 
here we might hope that this residual component is only a 
weak perturbation, the theory presented here is valid for 
any strength. This wide validity is of course very advan- 
tageous, however note that this approach is most useful 
in the uni-directional limit, i.e. when the residual terms 
are small in addition to being poorly phase matched jl] 1 
Here I only consider the effects of diffraction in any de- 
tail; other effects are not the specific subject of this work, 
and due to their lesser significance (here) are assumed to 
be incorporatd in the residual term Q. 



III. BEAMS, RODS, PIPES, AND DISKS 

The cylindrical geometry is perhaps the most likely to 
give useful results, as it covers not only the common case 
of a light beam of circular profile, but also propagation 
around the edge of a disk resonator. Here the coordinates 
are the axial z, the radial p and an angle 9. In the rest 
of this section, I choose each in turn as the direction 
for the underlying propagation, although by far the most 
common case is the axial "axi-symmetric" case relevant 
for the typical light beam. 



1 If the forward field has a wave vector fco evolving as exp(+«fcoz), 
the generated backward component will evolve as exp(— ikgz). 
This gives a very rapid relative oscillation exp(— 2ifcoz), which 
will quickly average to zero. 




FIG. 1: Longitudinal propagation in cylindrical coordinates: 
here the pulse propagates towards larger z, whilst transverse 
variation occurs in the radial p and angular (f> directions. 



To proceed we will need an expression for the Lapla- 
cian V 2 in cylindrical coordinates, which is just the usual 
expression 



1 



1 



V'E = -O p (pd p E) + —diE + diE, 



(5) 



where now E = E(z, p,c/)]uj). From this point we only 
need choose a primary propagation axis according to our 
interests, and proceed from there. In the following, I 
consider each possible choice in turn. 



A. Axial 

This axial case is suitable for the common case of free- 
space beam propagation, or that along a slowly changing 
rod or circular waveguide, such as a tapered optical fibre 
p"l[ . This is because we would expect angular variation 
and radial variation to be small and/or only slowly vary- 
ing. Note that axial propagation along the z coordinate 
has the nice feature that the propagation coordinate is 
translationally invariant along itself; i.e. we do not have 
to care where "z = 0" is. 

As indicated on fig. we choose the propagation direc- 
tion along the z axis of the cylindrical coordinates, with 
radial coordinate p and angular coordinate 9 to account 
for any transverse variation. This contrast with models 
(e.g. [y) which use the Cartesian x, y as transverse co- 
ordinates, although of course the use of radial transverse 
coordinates is far from unknown [refs]. We rewrite the 
wave equation (^|) to focus on z-propagation, and demote 
radial and angular effects to the status of residual terms, 
resulting in 

[d 2 + K 2 ] E = —Q - -d p pd p E - \d 2 E, (6) 
p p z v 

where the total wavevector is given by K 2 = n 2 (lu)uj 2 / c 2 . 
Factorizing gives us 

iQ 



d z E± = ±zKE ± ± — ± ^dppdp [E+ + E~] 



(7) 
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FIG. 2: Radial propagation in cylindrical coordinates: here 
the pulse propagates towards larger p, whilst transverse vari- 
ation occurs in the axial z and angular <j> directions. 



Here, in addition to the Q residual term, we have two 
additional coordinate-based residual terms. The first is 
that which gives radial diffraction, and the second an- 
gular diffraction. Both of these appear to have potential 
singularities at p = 0, but this is a coordinate effect - the 
singularity is not present in Cartesian coordinates. Thus, 
E will typically be smooth enough so that this will not 
cause pathological difficulties. 

Assuming both the radial and angular diffraction terms 
are small, we can decouple the E + and E~ fields as de- 
scribed and justified in more detail in H. For this to 
hold, we need all the residual terms on the RHS to be 
much smaller than the leading KE ± term, i.e. 



\dl (E+ 



|Q| < \2K 2 E ± \ 

)| < ^V^l 

I < \2K 2 p 2 E ± \ 



(8) 
(9) 
(10) 



These being sufficcntly well satisfied, we can approxi- 
mate eqn. (Q) to get the uni-directional wave equation 
for propagation in a beam or rod which is 



d z E ± = ±iKE ± ± 41 ± 



2K 2pK 



d p pd p E ± ± 



2 77± 



2p 2 K 



diE 



(11) 



B. Radial 

The radial case might be applied to the case where a 
wire or point source is radiating outwards into a cylinder 
or disk; or perhaps the reverse situation with converging 
fields. Alternatively, it might be useful when approach- 
ing the far-field, where part of an expanding wavefront 
enters some area of interest. Unlike the axial case where 
the absolute location z — was unimportant, here the 
coordinate centre at p = is fixed. 

As shown in fig. 0, we choose the propagation direc- 
tion along the p radial coordinate, with the axial z and 
angular 9 coordinates to account for any transverse vari- 
ation. We rewrite the wave equation ([|) to focus on p- 
propagation, and demote axial and angular effects to the 
status of residual terms, resulting in 



1 



-d p pd p + K 2 



E = Q - d 2 E - 4c)?£\ 



(12) 



where E = E(z, p, <j>; lo), and the total wavevector is given 
by K 2 = n 2 (uj)u} 2 / c 2 . Now since 

p-%pd p E = p- l d p [d p pE - E] , (13) 

then with F p = pE, we get 



[d 2 p+ K 2 ]Fp 



-pQ-pd 2 Fp-^d 2 F p + d p ^. (14) 



Factorizing gives us 
d p F± = ±iKFf ± 



ipQ 



2K 2K 



± 



ip 



91 K 



f: 



i F+ 



F: 



2K 



P 



(15) 

One feature of this is that the RHS has a residual term 
containing a d p derivative, which was generated when 
moving from the field E to the radially scaled version 
F p . We could therefore move this to the LHS now, but 
for simplicity I delay this adjustment until after the uni- 
directional approximation is made. 

Here, in addition to the generic Q residual term, we 
have three additional coordinate-based residual terms. 
The first is that which gives axial diffraction, and the 
second angular diffraction. The third arrives as a result 
of eqn. (n3J) , and acts as a radial drift. The second and 
third of these appear to have potential singularities at 
p = 0, but this is a coordinate effect - the singularity 
is not present in Cartesian coordinates. Thus, F will 
typically be smooth enough so that this will not cause 
pathological difficulties. 

Assuming both the axial and angular diffraction terms 
are small, we can decouple the E + and E~ fields as de- 
scribed and justified in more detail in For this to 
hold, we need all the residual terms on the RHS to be 
much smaller than the leading KE^ term, i.e. 

(16) 
(17) 
(18) 
(19) 



\ P Q\ < \2K 2 Ff\ 



\pdl 



f: 



< 1 | 



\2K 2 p 2 Ff\ 



\d pP - 1 {F++F p 



< |2K 2 F±| 



Of these, note in particular the last one, where we can 
see that we will need to be away from the origin for it to 
hold - as indeed might be expected on physical grounds. 
These being sufficently well satisfied, we can approximate 
eqn. @ to get the uni-directional wave equation for 
outward or inward radial propagation, which is 

tpQ 



dpFf = ±iKF± ± 



± ^d 2 l 
2K 2K 2 



± 



i F ± 
T 2K°" p ■ 
Now we can combine the two d p derivatives, 



(20) 
to get 



cU 1± 



2Kp 



±%KFf ± 



ipQ 



± 



?.p 



2K 2K 



d 2 F± ± 



(21) 
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FIG. 3: Angular propagation in cylindrical coordinates: here 
the pulse propagates around in positive (/>, whilst transverse 
variation occurs in the axial z and radial p directions. 



C. Angular 

Angular propagation is relevant where the light is prop- 
agating around some kind of circular waveguide, such 
as in a whispering-gallery (disk) waveguide, although it 
could also be applied to a helical waveguide. The restric- 
tion to waveguides results from the fact that without 
some confining structure, light will travel in a straight 
line, and so would only only be nearly angular for a brief 
interval at closest approach to the coordinate origin. The 
angular case has the nice feature that the propagation 
coordinate (?) is translationally invariant along itself 
(around the origin); i.e. we do not have to care where 
"0 = 0" is. 

As shown in fig. |3|, here we choose the propagation 
direction around the (f> angular coordinates with the axial 
z and radial 9 coordinates to account for any transverse 
variation. We rewrite the wave equation (^) to focus on 
(^-propagation, and demote axial and radial effects to the 
status of residual terms, resulting in 



Id 2 



E = Q- d 2 E - -d p pd p E, (22) 



with E = E{z, p, 4>; uj), and where the total wavevector 



is given by K 2 
get 



! (w)w 2 /c 2 . Then with F^ = p 2 E, we 



1 



[dl + K 2 ] F = p 2 Q - p 2 d 2 F p - -d p pd p F p . (23) 



Factorizing gives us 



ip 2 Q 



% P 2 a2 



B*Ff = ± ip KFf ±^±£=8 



± ^d pP dp 



2K 2K 

Ft + f: 



F+ + F- 



(24) 



Here, in addition to the Q residual term, we have two 
additional coordinate-based residual terms. The first is 
that which gives axial diffraction, and the second radial 
diffraction; although the second may be split into an al- 
ternate diffraction along with a radial drift by using eqn. 



Assuming both the axial and radial diffraction terms 
are small, we can decouple the E + and E~ fields as de- 
scribed and justified in more detail in 0|. For this to 
hold, we need all the residual terms on the RHS to be 
much smaller than the leading KE^ term, i.e. 



IpQ| « 



F 



F: 



d P pdp F+ + F, 



2K 2 Ft 



2K 2 F± 



2K 2 F± 



(25) 
(26) 
(27) 



These being sufficently well satisfied, we can approximate 
eqn. ( p4| ) to get the uni-directional wave equation for 
angular propagation, which is 



= ±tpKF± ± 
ip 



ip 2 Q 
2K 



± l -^d 2 F ± 
± 2K° zt * 



± 7^d pP dpF i 



2K 



(28) 



Since to maintain uni-directionality during this kind 
of angular propagation, our wave must be somehow con- 
fined in a ring shaped waveguide, most likely the radial 
terms included above will not be relevant - any radial 
diffraction will have been already balanced by the radi- 
ally confining waveguide structure, and the radial wave 
profile will match some guided mode. In this case, we 
can use 



d^Ff = ±i P KF± ± 



^±^9 2 F 
2K 2K z < 



± 



(29) 



to propagate light pulses around a thick disk, ring, or 
pipe, whilst still allowing for axial diffraction. 



IV. CONCLUSION 

Here I have derived bi-directional factorizations of the 
Helmhotz wave equation in the cylindrical geometry, fo- 
cussing on each possible choice of propagation direction 
in turn. These then allow approximate uni-directional 
forms, based on a generalized notion of paraxiality; and 
it is these which are likely to be most useful. These re- 
sults are done in the same style as, and are intended 
to complement existing calculations done using cartesian 
coordinates [Q. 

One could certainly also image following this same 
procedure using other orthogonal coordinate systems 2 , 
notably spherical-polars or parabolic coordinates. You 
could also use the approach to model diffraction of ray- 
like light beams in a conformal cloak 16 or similar; pro- 
viding sufficient physical motivation exists and the a uni- 
directional approximation to the resulting wave propaga- 
tion equations is achievable. 
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Appendix: Factorizing 

Here is a quick derivation of the factorization process; 
the z-derivative has been converted to ik, 1 — n 2 uj 2 /c 2 , 
and the unspecified residual term is denoted Q. 



[-k 2 +p 2 ]E=-Q 



E 



k 2 ~^ Q 



(k - P) (k + /?) 



2/3 



Q. 



(30) 
(31) 

(32) 



_k + (3 k-f3_ 
Now (k — /3) _1 is a forward-like propagator for the field, 
and (k + /3)" 1 a backward- like propagator. Hence write 
E = E + + E~ , and split the two sides up 



2/3 



k + /3 k~/3 



Q 



E ± = 



[k T 13} E ± 
ikE ± 



— — !— Q 

2/3 

±^E± ± 



(33) 
(34) 
(35) 
(36) 



and reverting to z derivatives gives us the final form 



dzE* = ±i(3E ± ± —Q. 



(37) 
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